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We introduce a one- dimensional sandpile model which in- 
corporates particle inertia. The inertial dynamics are gov- 
erned by a new parameter which, as it passes though a thresh- 
old value, alters the toppling dynamics in such a way that 
the system no longer evolves to a self-organized critical state. 
A range of mean-field theories based on a kinetic equation 
approach is presented which confirm the numerical findings. 
We conclude by considering the physical applications of this 
model, particularly with reference to recent experimental re- 
sults. 

PACS numbers: 64.60.Ht, 02.50.-fs, 05.40.+j, 05.60,+w 



I. INTRODUCTION 

The concept of self-organized criticality (SOC) was in- 
troduced by Bak, Tang and Wiesenfcld (j]] as a pos- 
sible explanation for the common occurrence of scale- 
invariance in nature. To demonstrate this behaviour, 
they introduced the sandpile model, a driven dissipative 
cellular automata whose dynamics are defined by local 
interactions. Despite the short-range dynamics, the sys- 
tem organizes itself into a non-equilibrium critical state 
with no finite correlation length and hence no character- 
istic length scale. A feature common to all SOC systems 
is that the sizes of fluctuations follow power-law distribu- 
tions, a direct consequence of the scale-invariance. How- 
ever, comparisons with real systems have met with only 
partial success. Power-laws were observed in a granular 
mixture when avalanches were initiated by watering the 
pile H ,but not in a pile of glass beads that was gradually 
tilted Adding grains individually to a conical sand- 
pile only showed power-laws for sufficiently small piles [Q . 
Recently, Frette et. al. || performed experiments on one- 
dimensional piles of rice and found power-law behaviour 
only for grains with sufficiently large aspect ratio, fn 
light of these experiments, we believe it would be infor- 
mative to construct a sandpile model with greater phys- 
ical applicability, whilst hopefully retaining some of the 
interesting dynamics. 

One important ingredient missing from the basic model 
is inertia. To the best of our knowledge, there have 
only been two attempts to construct an inertial sandpile 
model. Prado and Olami Q chose to associate moving 
particles with a decrease in the local stability, and found 
SOC behaviour only for small systems, in both one and 
two dimensions, thus giving an explanation for the re- 
sults in H . Krug, Socolar and Grinstein j7| gave a single 
measure of momentum to the entire cluster of moving 



particles in a one-dimensional system. They found that 
their inertia parameter needed to be fine-tuned to zero 
for the system to become SOC. 

In this paper, we consider a new way of incorporating 
inertia into the sandpile model. A full description of the 
model is given in the following section, but briefly, we 
suppose that moving particles only come to rest on those 
sites whose slopes are not too steep, where the definition 
of 'too steep' is controlled by a new parameter. Our 
main result is the crossover between SOC and non-SOC 
behaviour as this parameter passes through a threshold 
value. The existence and value of this transition point has 
been confirmed by extensive numerical analysis, mean- 
field theories and qualitative reasoning. 

This paper is arranged as follows. The new model is 
defined and some of its i mm ediate consequences explored 
in section O. In section III , the results of numerical sim- 
ulations are given for the simplest non-trivial sandpiles. 
These results are confirmed by the mean-field analysis 
given in section IV , where a rate equation approach has 
been adopted. Finally, in section [v|, we explain the na- 
ture of the transition threshold and discuss applying the 
model to real physical systems. 



II. THE MODEL 

A one-dimensional sandpile is defined by a set of in- 
teger heights hi, i — 1 ... L, or equivalently by the local 
slopes Zi — hi — hi + \. The right-hand boundary is taken 
to be open, Zl = hi,, whereas the left-hand boundary 
is treated as closed, z = 0. Particles are added sequen- 
tially to randomly chosen sites, increasing their height by 
one unit. In the limited local sandpile (LLS) model [||, 
a site i becomes unstable when Zi becomes greater than 
the critical slope parameter z c . Any such unstable site 
will topple, z c particles leave site i and move to site i + 1 
(or leave the system if i = L), so 

Z%—\ ► Zi—l -f- Z c , 

Zi ► Zi 2z c , 
Zi+i -> Zi+i + z c , (1) 

with equivalent rules for toppling at boundaries. It is 
now possible for Zi—i and/or Zj+x to become unstable 
and topple, and an avalanche will begin. We call the 
series of the initial topplings at sites i — 1, i — 2, i — 3 
...the back- avalanche. A back-avalanche propagates to 
the first site j < i with Zj < 0, where such sites are called 
troughs. In terms of the slopes before particle addition 
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to site i, Zi-i — 1 will also stop a back-avalanche. Some 
previous literature || refers to such instances as slide 
events. Topplings occur on a timescale much smaller and 
separated from that of particle addition, in that no more 
particles are added to the pile until the avalanche is over 
and every site has become stable (z.; < z c Vi). 

As this model stands, the inertia of the toppling par- 
ticles has been ignored, any momentum in the topple is 
assumed to be dissipated instantaneously from the sand- 
pile as soon as the particles have moved. The previ- 
ous attempts to incorporate inertia into the system [^]|7]] 
used quantities that evolved throughout the avalanche 
and hence introduced a form of memory into the sys- 
tem. In this paper, we consider a new set of rules for 
toppling that requires just one extra, time-independent 
quantity, but nonetheless intuitively mimics inertia. We 
introduce the minimum- slope inertial limited local sand- 
pile (MILLS) model, which has a second critical slope 
parameter, the inertial threshold z in . Now when a site i 
topples, the particles will be deposited on the first site 
j > i obeying 

Zj <C Zi n , (2) 

or (if no such site exists) leave the system. This calcula- 
tion is performed right-to-left for each unstable pile and 
individually for each of the z c particles involved in every 
toppling, recalculating the slopes as a particle comes to 
rest. An example is given in fig.|l] for z c = 2 and Zi n = 3. 
Note that we have implicitly introduced a timeframe for 
inertial effects which is much smaller than that of top- 
pling events. This separation of timescales means that 
the time taken for an unstable site to begin to topple is 
much longer than that for the ensuing particle motion 
itself. 

As a realistic physical system, this model does have 
some obvious drawbacks, the most notable being that a 
particle will be stopped just as easily if it has just top- 
pled, or if it has rolled down a large region of slopes 
> Zi n . We could of course make z,- n dependent on each 
particle's previous motion, but it is the lack of inher- 
ent memory that makes much of the subsequent analysis 
possible. Furthermore, with this implementation of iner- 
tia the sandpile takes a microscopically smooth, realistic 
shape which further justifies its study. 

We can immediately make some general observations 
about the MILLS model. Since the particles toppling 
from an unstable site all move at least one step, the 
slope of the site to the immediate left will always increase 
by z c independent of the value of zi n . Thus, the back- 
avalanches in this model are identical to those in the stan- 
dard model, and consequently troughs are still important 
as bounds of the left-hand edge of an avalanche. How- 
ever, troughs no longer necessarily bound the right-hand 
edge of an avalanche as they did in the LLS model. In- 
deed, no expression involving a single site can now serve 
as a general right-hand bound, so the previous analysis 
of the standard model based on troughs |7|,[l0) cannot be 
extended to this model. 



For z c = 1, the sandpile soon reaches a trivial steady 
state with = 1 Vi independently of z in . Inertial effects 
are now indistinguishable from topplings, added parti- 
cles move downslope by either mechanism and leave the 
system in an unaltered state. Varying Zi n will change 
the time taken for the particle to reach the right-hand 
boundary, but only for z c > 2 can Zi n have any influence 
on the nature of the critical state. 

If Zi n = — oo, inertial effects will dominate and all top- 
pling particles will immediately leave the system, as if the 
sandpile has become 'infinitely slippery'. This is equiv- 
alent to the case r = for the model in [j. Since sites 
with negative slope are rare, Zi n < will have a similar 
effect to Zi n ~ —oo. For Zi n < z c , no particle will come 
to rest on a site with critical slope and so no sites to the 
right of the initial topple can become unstable. As a re- 
sult, there will be no forward toppling and the avalanche 
will consist solely of the back-avalanche. 

It may appear that implementing inertia in this man- 
ner could allow for a large cluster of sliding particles to 
all come to rest simultaneously on the same site, say a 
site i with z% large and negative. It is easy to prove, 
however, that at most z c particles can stop on a site dur- 
ing any one step in the toppling process. To see this, 
first realise that a back-avalanche consists of one topple 
per avalanche time step, so to get more than z c particles 
moving at once we must have Zi n > z c . Furthermore, a 
site k < i must topple and subsequently allow at least 
one particle to slide through it, ie. 

Zk Z c > Zin, (3) 

so for this sequence of events to occur requires z^ > 2z c . 
Although site k + 1 toppling could contribute to such 
a large slope, having Zk+i > z c on the previous step 
would simultaneously require Zk < 0, except at the left- 
hand boundary of the avalanche, as closer analysis soon 
reveals. It should now be apparent that for any site to 
gain more than z c particles during any one avalanche 
step, there must already exist a site in the system to 
which this has already happened. Since particle addition 
cannot create such a site, we must conclude that they 
never occur. Intuitively, this result implies that inertial 
effects serve only to broaden the shape of the avalanche. 

If Zi n = oo, inertial effects will be so weak that toppling 
particles will never move by more than one site, and we 
just have the LLS model. In fact, from the result given 
in the previous paragraph, it should be apparent that the 
maximum slope that a site can reach is 2z c , so if Zi n > 
2z c , inertial effects become redundant and the MILLS 
model reduces to the standard model. 
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might relate to actual systems. Suppose the particles all 
have an aspect ratio a. A vertical stack consisting of 
rounded particles will clearly be much less stable than a 
stack of flat particles, so z c will decrease as a decreases. 
The parameter z;„ should also decrease with a, since 
rounder particles will roll more easily and each site will 
be narrower and hence easier to traverse. Indeed, varying 
a should to change faster than z c does, so we 

can see that the ratio ' z c will increase with a. More 
precise analysis in this manner is possible but any such 
accuracy is lost within the artificial framework common 
to all sandpile models. 

In summary, the MILLS model differs from the LLS 
model by the extra parameter Zi n - For Zi n > 2z c we 
just get the LLS model, for < z c the toppling process 
consists of just the back-avalanche, and for < the 
majority of sliding particles leave the system immediately 
after toppling. 



w 

WW, 






w 





































(b) 



(c) 



FIG. 1. An example of the toppling process in the MILLS 
model, for z c = 2 and Zi n = 2. (a) Particle added to a site 
with critical slope. The site becomes unstable and topples, 
(b) The first toppling particle sticks to the adjacent site, giv- 
ing it a local slope of 3. The site marked with an asterisk is 
still in the process of toppling, (c) The second particle slides 
throgh the site with slope 3 to stop at the next site. All sites 
in this region are now stable. 

In (|] it was suggested that it was the shape of the 
particles that determined the nature of the critical pile. 
Although any sandpile model is over-simplified compared 
to real piles of any substance, for instance in failing to 
allow for the particles to overlap and for ignoring any 
variation in particle dimensions, it is still instructive to 
consider qualitatively how the two parameters z c and Zi n 



III. RESULTS 

Numerical analysis of finite one-dimensional sandpile 
models is difficult because the convergence to the asymp- 
totic SOC regime is very slow. Indeed, we are not 
aware of any simulations sufficiently large to demonstrate 
the power-law behaviour expected in the thermodynamic 
limit, so an alternative test for self-organised criticality 
must be employed. If the critical state is governed by 
exponentially decaying correlations, then increasing the 
system size far beyond the correlation length could not 
alter the avalanche size. However, if instead the correla- 
tions decay as power-laws, then the critical state is scale- 
invariant and so the largest avalanche sizes will always 
vary with the system size. Two examples of dropsize fre- 
quency distributions are given in fig.|^, where the dropsize 
is defined as the number of particles to leave the pile as 
the result of a single particle addition. The distribution 
broadens for Zi n — 2z c when the system size is doubled, 
which is clearly not the case for z^„ = 0, thus the system 
is not SOC for Zi n = 0. It is obviously important to find 
the value of Zi n between these two extremes at which the 
transition between SOC and non-SOC behaviour occurs. 

The distribution of dropsizes is not always a useful 
measure of avalanche size, since for Zj„ = z c the max- 
imum dropsize is just z c + 1, as the following analysis 
demonstrates. Suppose an avalanche is started by the 
addition of a particle onto a site i with critical slope, so 
that before toppling we have z% = z c + l and z;_i < z c — 1. 
When site i topples, Z{ will decrease by at least z c , pos- 
sibly more if any particles stick to i + 1, so now zi < 1 
and Zi-\ < 2z c — l. Presuming now that i — 1 topples, all 
but possibly one particle will stick to i, giving < 0. 
Even if i — 2 now topples, no particles can slide past 
i — 1 and so now z^_2 < 0. From now on this is all 
that happens, toppled particles neither slide nor cause 
any further topples. The avalanche will eventually end 
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with at most z c + 1 particles from the leading edge of 
the avalanche moving beyond site i and potentially con- 
tributing to the dropsize. More qualitatively, toppling is 
minimal for z;„ < z c and the bulk of the sliding is limited 
to within the avalanche for Zi„ > z c , so for Zi n — z c we 
should expect particle transport to be low. Less com- 
monly considered measures of avalanche size, such as the 
total number of topples involved, must be used in this 
case. 
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FIG. 3. Dropsize distributions for z% n = 1 and z c = 2, 
with varying system size. Asterisks correspond to L=128, 
open circles to L=256, crosses to L=512 and filled circles to 
L=1024. 



Log( Dropsize = Number of particles to leave pile ) 

FIG. 2. Distribution of dropsizes for Zi n — 4 and 0, and for 
two different system sizes, with z c = 2. Circles correspond to 
L = 512, crosses to L — 1024. 

For Zi n — 1 and z c — 2, the dropsize distributions 
broaden with the system size, as shown in fig.[| There 
is perhaps some hint of a lessening in the broadening for 
the larger systems. However, long processing times have 
resulted in statistically noisy data which it is difficult to 
analyse. An alternative does exist, since when z,;„ < z c 
the avalanche size is bounded above by the distance be- 
tween troughs, so a system with a finite concentration of 
troughs must have a finite correlation length and hence 
cannot display SOC. The trough density, being just a 
single value, is much less susceptible to noise and can 
be measured reliably for larger L. The distribution of 
trough densities for z c — 2 and for L up to 2048 is given 
in fig^. That Zi n = 1 is not SOC is now evident as 
the trough density tends to a finite value, although this 
value is small, corresponding to a large correlation length, 
which explains the observed broadening of the dropsize 
distribution in small systems. For z;„ = z c = 2, the 
trough density approaches zero with system size faster 
than for any other value of Zi n , corresponding to a singu- 
lar correlation length and SOC behaviour. These results 
are borne out by the mean-field analysis in section IV. 
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Another useful measure to consider is the average slope 
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where tan _1 (S') is the angle of repose. The variation of 
S with L and Zi n is given in fig.^[ For Zj„ < 2 the lines 
converge on values in the range 1 < S < 1.2. It is already 
known |ll]] that for Zi n ~ 4 the slope for L — > oo is S = | . 
Comparing this line to those for other Zi n in figj5[ it could 
be judged that the curves for Zi n — 2 and Zi n — 3 are 
both tending to S = 2 in a similar manner. Since conver- 
gence is slow, however, verifying these asymptotic limits 
is impossible from this data alone. Taking logarithmic 
plots to project the lines further, as was done success- 
fully for the standard model jy], fails here as it gives 
reasonable straight line fits for a range of limiting slopes. 
Instead, a rough prediction of the asymptotic slope for 
Zi n — 2 is given by the following qualitative argument. 

For Zin < z c particles leaving the pile cannot cause 
the rightmost site i = L to topple, although if zl < Zi n 
they will stick there and increase Zl- Once zl is in the 
range Zi n < zl < z c , all particles will slide over L, so zl 
is fixed until the end of the current avalanche. The same 
thing will then happen for site L — 1, then L — 2, and 
so on. Thus, the effect of any avalanche reaching the 
right-hand boundary is to leave behind a cluster of sites 
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with slopes in the interval [zi n ,z c ]. Furthermore, the 
slopes in this region are stable under subsequent particle 
addition on any site to its left. Similar clusters should 
grow throughout the pile for similar reasons, although 
their right-hand edges will not be fixed. This has special 
significance for Zj„ = z c , since the clusters of sites will all 
have slope z c and S is expected to tend to a value at least 
close to z c . Given that this system becomes dominated 
by sites of critical slope, the z c + 1 particles that slide 
from any avalanche should all leave the pile, thus the 
modal dropsize should also be the maximum one, z c +l. 
This has been confirmed by the numerical studies. 

i i i i i I i i i i I i i i I i i i I 



is varied. In this section we present a range of models for 
Zi n < 4 with z c = 2, comparing their predictions to the 
numerical findings. All these models are cast in terms 
of the slope distribution {S n (t)} ; where each S n (t) is de- 
fined to be the proportion of sites in the system at time 
t with slope n. For convenience, we also define 

o 

S T (t) = ]T Si(t), (5) 

i— — oo 

which is the trough density. As these slopes represent 
the sandpile in its stable state, the maximum value any 
Zi can take is two and 
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FIG. 4. Trough densities for z c = 2 against system size. 
Circles refer to Zi n = 4, squares to Zi n = 3, diamonds to 
Zi n = 2, plus-signs to Zi n = 1, triangles to Zi n = and crosses 

tO Zi n = —OO. 



IV. MEAN-FIELD ANALYSIS 

Christensen and Olami successfully formulated a 
random nearest neighbour mean-field model for critical 
height sandpiles. Directly extending their procedure to 
critical slope models has proved to be intractable, but the 
principle of taking random neighbours in an infinite sys- 
tem can still be applied to the MILLS model with some 
success. Utilising this approach is non-trivial, however, 
as the dynamics of the system change significantly as Zi n 



Sn 



S 1 + S 2 = 1. 



(6) 



The simplest case to analyse is Zi n = — oo, where any 
avalanche propagates back to the first trough and all top- 
pled particles are removed from the system. Particle ad- 
dition alone has the effect of increasing the slope of a 
randomly chosen site from n to n + 1 , 



Sn > S n — , 
S n +1 — > SVi+1 + £ 5 



(7) 
(8) 



and in a similar manner decreasing the slope of another 
randomly chosen site by one. Troughs are assumed to be 
primarily sites of slope zero, so increasing the slope of a 
trough will give a site of slope one, whereas the effects of 
decreasing the slopes of troughs are ignored. 

When a particle is added to a site of slope two, (@-||) 
must be replaced by equivalent rules for the whole of the 
subsequent sequence of topples. There are two cases to 
consider - if the site to the immediate left (here chosen at 
random) has slope one, a slide event occurs, removing two 
particles from the system but leaving {S n (t)} unchanged. 
Otherwise, we get a full avalanche which only alters the 
slopes at its leading edge (where the particle was added) 
and its trailing edge (at the trough that halted the top- 
pling). The net result of both these edge effects is the 
lose a trough and a site of slope two, and to gain two sites 
of slope one. It is now straightforward to write down rate 
equations for {S n (t)}, with the timescale normalised to 
L particle additions per unit t, 



dS 2 
dt 

dSi 
dt 

dSx 

~dT 



-2S 2 + 5i + S x Sa, 
—2S\ + St + 35*2 — SiS 2 , 
—St + Si — S 2 . 



(9) 
(10) 

(11) 
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These equations can be justified directly - for instance, 
the first term on the right hand side of (]L0|) corresponds 
to particle addition onto or adjacent to a site of slope one, 
the second term is for a particle added to a trough, the 
third term accounts for particle addition adjacent to or 
onto a site of slope two, where the latter instance invokes 
an avalanche, and finally the fourth term subtracts out 
trivial slide events. In practice, however, it is simpler to 
derive these equations by writing down all nine possible 
combinations of pairs of sites, ie. (TT) (Tl) (T2) (IT) 
(11) (12) (2T) (21) and (22), then computing directly 
the effects of particle addition to the right hand site and 
weighting them accordingly. 

The equations (|9|-|ll|) converge on the stable solution 
(St, Si, S 2 ) — (|, 3, ¥J, in good agreement with the nu- 
merical values (0.18,0.48,0.34). Since the influence of 
sites with strictly negative slope have been ignored, this 
mean- field theory actually covers the range of all Zj„ < 0, 
and indeed the agreement still holds up to Zi n = 0. For 
zin = 1 the avalanching particles can no longer be ig- 
nored and new terms need to be added to account for 
where they come to rest. The expected number of par- 
ticles brought into motion during an avalanche depends 



on the state of the site to the immediate left of the site 
where the first particle was added. If this site has slope 
less than two then just the one site topples. Otherwise, 
the number of sites that topple is half the average dis- 
tance between troughs. In this random neighbour model, 
the probability that two nearest troughs are a distance 
r apart is (1 — St) t St- Consequently, half the mean 
distance between troughs is (1 — St)/2St and the total 
number of sites that topple during an avalanche is on 
average S T + Si + S 2 [(1 - S T )/2S T ]- Each toppled site 
releases two particles, so we will need a factor of two. 
Furthermore, since avalanches are initiated by particle 
addition to sites of slope two, we will also need a factor 
of S < 2. Thus, the expected number of toppling particles 
per unit time is 



A(t) = 5 2 ( 2(1 - S 2 ) + ^-^S 2 

dt 



(12) 



Each of these X(t) particles will stop on a trough, creating 
a site with slope one. This will also cause the site to the 
immediate left to decrease its slope by one. Moreover, we 
know that the particle passed over this site, so it cannot 
be a trough. Therefore it must have slope one or two 
with probability Si/ (I — St) or 52/(1 — St), respectively. 
Since we are ignoring any correlations between where a 
particle starts moving and where it stops, all we need to 
do is to add extra terms to the right hand sides of (p|-|ll|) 
to account for the changes in the slopes caused by the 
X(t) particles coming to rest, 



— r— — —25 2 + 5i + 5i5 2 — A 



d5 2 
dt 
dSi 
~dT 

cLSt 
~dT 



I -S T J ' 



c 



-2N, - .S, f- :-!.S' : , - ,s'| ,s'_> + A ( I + f 2 ) , ( 



1-S T 



Si — S 2 + A 



Si 



1-St 



c 



These equations now converge to (5t,5i,5 2 ) ~ 
(0.09, 0.65, 0.26), whereas numerically the slope distribu- 
tion is (0.07, 0.68, 0.25). Note that the trough density 5 T 
is small but non-zero, corresponding to a large but finite 
correlation length. Thus, as noted in the previous sec- 
tion, large system sizes are required to demonstrate that 
this system is not SOC. 

For Zi n — 2, the unique nature of the avalanche re- 
quires a new set of rate equations, which can be derived 
in a similar manner to before by considering the conse- 
quences of particle addition to a pair of sites with arbi- 
trary slopes, 



§ + * ♦/•(lis 



1 



dSi 
~dT 

dST 
~dT 



-25? + 5 T (35 2 - Si + S T ) + n 1 



St — Si 

l-S 2 



-St + Si — S 2 Si — [i 



St 



I-S2 
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where fi(t) = S 2 (2 + S 2 ) is the expected number of 
particles that emerge to the right of the particle addi- 
tion. Although these equations admit the steady solution 
(S T , S u S 2 ) -> (0, 0, 1) with Si ~ 1 - £2 and S T ~ 0, cu- 
bic corrections are required for it to be stable. Numeri- 
cally, the slopes are (0.01, 0.08, 0.91) for L = 2048 but, as 
mentioned in section III , the asymptotic slope appears to 
be S = 2 which demands the slope distribution (0, 0, 1). 

For Zi n = 3, all the sites involved in an avalanche 
resolve in a highly non-trivial way, so calculating rates 
between the states before and after an avalanche be- 
comes impossible. Instead, we introduce the following 
local slope argument to predict the slope distribution in 
the critical state. Within an avalanche, a toppled site 
deposits a pair of particles to its right, which, after any 
inertial motion, settle and possibly cause further top- 
pling. Keeping our random neighbour approach, the ex- 
pected number of further topples caused takes the same 
constant value E t throughout the avalanche. If E t < 1 
then avalanches would die out exponentially quickly, par- 
ticles will not get transported from the pile and so the 
sandpile would build up, increasing the average slope and 
hence also E t . If Et > 1, avalanches would grow expo- 
nentially and large numbers of particles would leave the 
pile, decreasing E t . Although transport from the pile 
is not actually catered for by this model, we infer that 
E t = 1 in the critical state, so finding an expression for 
E t in terms of {S n (t — 00)} gives an additional constraint 
which might help to fix the slope distribution. 

This local slope argument has several deficiencies, 
apart from the obvious drawbacks of the random neigh- 
bour approach. Back-avalanches have been ignored, so 
this method will certainly break down for z,- n < 2. The 
effects of single particle addition has not been accounted 
for, which will further exasperate the situation when S2 
is small. Although this gives a limited range of applica- 
tions for this model, it has the advantage of being con- 
ceptually and mathematically simple. For instance, in 
the standard model z.- ln = 4, adding two particles to any 
site with a slope greater than zero will cause one further 
topple, so Et — 1 — St- Requiring that this equals one 
gives St = 0, as expected. 

For Zi n — 3, consider the effect of adding the particles 
to two adjacent sites, i and i + 1. If Zi < 2, then both 
particles would stick to i, giving the site a slope Zi + 2, 
which will topple if Zi = 1. However, if Zi — 2 then 
the last particle added will slide onto the adjacent site, 
resulting in Zj+i — » Zj+i + 1 and Zi remaining unchanged. 
Thus, the expected number of further topples is 

E t =Si + Si (19) 

Setting this equal to one and employing (^), we get an 
equation for the trough density, 

S T = S 2 (S 2 - 1). (20) 

By definition, all the Si must lie in the range [0,1]. How- 
ever, (p0) is strictly negative for < S2 < 1, so for 



consistency we must have S2 = or 1. Avalanches could 
never begin in a system with S 2 = 0, so we conclude that 
(St, Si, S 2 ) — (0, 0, 1), in agreement with the numerical 
work. 



V. CONCLUSIONS AND DISCUSSION 

It has been demonstrated that the MILLS model with 
z c — 2 displays self-organized criticality for systems with 



> z c , but SOC is lost when Zi„ < 



Zr. 



The SOC 



state for z,„ > z c was identified by the broadening of the 
dropsizc distribution with the system size, whereas for 
Zin < z c the existence of a finite correlation length was 
found both numerically and from the mean- field theories. 
Preliminary numerical results for z c up to 4 demonstrate 
that the transition point Zi n = z c holds more generally. 
We postulate that for all z c > 1 the model is SOC only 
when z^ > z c , and now argue why this should be so. 

An essential requirement of any SOC system is that 
it is governed by local driving forces For Zj„ > z c , 
inertial effects only come into play inside an avalanche, 
since it is only within an avalanche that sites with slopes 
greater than z c can occur. Thus, the couplings re- 
main short-range and SOC behaviour is preserved. For 
Zin < z c , however, particles can slide far beyond the 
avalanche, so interactions become non-local and the sys- 
tem ceases to be SOC. We might navely expect that a 
system with Zj„ = z c will also fail to be SOC, especially 
since the system becomes dominated by sites of slope z ci 
which should allow particles to slide through large dis- 
tances. However, the exact analysis in section [II proved 
that only a tiny proportion of particles move beyond the 
avalanche. In fact, the vast majority move by just one 
site, which is certainly a local driving force, so we can 
conclude that a system with z.i n — z c is SOC. 

In the experiments carried out on piles of rice ||, 
power-law behaviour was only observed for rice with suf- 
ficiently large aspect ratio. In section [n] it was argued 
that Zinjz c should increase with the particle aspect ra- 
tio. Thus, if the aspect ratio is increased such that Zi n 
becomes greater than z c , the MILLS model gives a possi- 
ble explanation for the appearance of power-laws - the 
system has crossed over into the SOC regime. How- 
ever, many non-SOC systems also exhibit power-law be- 
haviour, and an alternative model has recently been pro- 
posed by Newman and Sneppen |0|| . Here, they claim 
that an SOC approach is invalid when the power-laws 
have characteristic exponents in the vicinity of 2, as 
was observed for the rice piles. Their model is instead 
driven by external noise, but still predicts a crossover 
from power-law behaviour. 

Regardless of its experimental applications, the MILLS 
model represents an interesting addition to the range of 
sandpile systems already studied, in that it provides an 
inertial sandpile model which is SOC for a broad band of 
parameter space - that is, without the need for fine tun- 
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ing - and independent of the system size. Extending this 
model to higher dimensions or to critical height sand- 
piles should in principle be straightforward, but since 
our original motivation was to try and account for the 
results from the rice pile experiments, we chose to base 
our model on that system. 
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